#R version 4.4.3 (2025-02-28 ucrt)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 11 x64 

library(dplyr) #1.1.4
library(ggplot2) #3.5.1
library(readr) #2.1.5
library(haven) #2.5.4
library(estimatr) #1.0.6
library(broom.mixed) #0.2.9.6
library(lme4) #1.1-37
library(countrycode) #1.6.1
library(forcats) #1.0.0
library(stargazer) #5.2.3
library(tibble) #3.2.1
library(modelsummary) #2.3.0

`%notin%` <- function(x,y) !(x %in% y) 

setwd("C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/cps/replication")

wth <- read_csv("data/wth_panel.csv") |>  
  select(countryname, year, v2xps_party, v2x_polyarchy, wth_duration, wth_regime, wth_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, wth_pduration, wth_number, fsu, ussr_sat, regime_elecs, av_dm) |>  
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown")) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(wth_prior %in% c("Multi-party", "Military", "One-party", "Monarchy") & is.nan(.) | wth_prior %in% c("Multi-party", "Military", "One-party", "Monarchy") & is.na(.), 0, .)) |>  
  filter(wth_prior %in% c("One-party", "Military", "Multi-party") & wth_regime == "Democracy") |>  
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) |>  
  mutate(wth_prior = factor(wth_prior, levels = c("One-party", "Multi-party", "Military"))) |>  
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) |> 
  group_by(ccode) |> 
  mutate(ldv = lag(v2xps_party), reactive = case_when(ruling == 1 ~ 0, ruling == 0 ~ 1)) 

#Figure 1.1 of Appendix C----
mp <- lm_robust(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + regime_elecs + av_dm, data = wth, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "Pooled") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "wth_priorOne-party", "wth_priorMulti-party", "wth_priorMilitary"))

df1 <- filter(wth, elections == "None")

m0 <- lm_robust(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "None") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "wth_priorOne-party", "wth_priorMulti-party", "wth_priorMilitary"))

df1 <- filter(wth, elections == "First")

m1 <- lm_robust(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "First") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "wth_priorOne-party", "wth_priorMulti-party", "wth_priorMilitary"))

df1 <- filter(wth, elections == "Second")

m2 <- lm_robust(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "Second") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "wth_priorOne-party", "wth_priorMulti-party", "wth_priorMilitary"))

df1 <- filter(wth, elections == "Third-plus")

m3 <- lm_robust(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + av_dm, data = df1, se_type = "HC1") |>
  tidy() |>
  mutate(Sample = "Third-plus") |>
  filter(term %in% c("pi1", "prior_pi", "asp_dummy", "wth_priorOne-party", "wth_priorMulti-party", "wth_priorMilitary"))

#fig1 left hand side
rbind(m0, m1, m2, m3, mp) |>
  mutate(Sample = factor(Sample, levels = c("None", "First", "Second", "Third-plus", "Pooled"))) |>
  mutate(var = recode(term, "pi1" = "Incumbent PI", "prior_pi" = "Prior PI", "asp_dummy" = "ASP")) |>
  filter(var %in% c("Incumbent PI", "Prior PI", "ASP")) |>
  ggplot(aes(x = Sample, y = estimate, ymin = conf.low, ymax = conf.high, color = var, shape = var)) +
  geom_pointrange(position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(y = estimate - 100), linetype = "solid", color = NA) +  # Setting color to NA to remove legend
  geom_point(aes(y = estimate - 100), shape = NA, color = NA) +  # Setting color and shape to NA to remove legend
  theme_bw() +
  xlab("") +
  ylab("") +
  ylim(-0.55, 0.75) +
  scale_color_grey(name = "Variable", start = 0, end = 0.8) +
  scale_shape_manual("Variable", values = c(15, 16, 17)) + 
  theme(legend.position = "bottom")  +
  ggtitle("OLS") +
  theme(plot.title = element_text(hjust = 0.5))

#ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_fig1_1.jpg", dpi = 500, width = 5, height = 5)

#Figure 1.2 of Appendix C-----
wth8 <- wth |> select(v2xps_party, wth_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, wth_pduration, wth_number, fsu, ussr_sat, regime_elecs, av_dm, elections, ccode) |> na.omit()

m1 <- lmer(v2xps_party ~ wth_prior + pi1*elections + prior_pi + asp_dummy + (1 | ccode), data = wth8, REML=FALSE)
m2 <- lmer(v2xps_party ~ wth_prior + pi1*relevel(elections, ref = "First") + prior_pi + asp_dummy + (1 | ccode), data = wth8, REML=FALSE)
m3 <- lmer(v2xps_party ~ wth_prior + pi1*relevel(elections, ref = "Second") + prior_pi + asp_dummy + (1 | ccode), data = wth8, REML=FALSE)
m4 <- lmer(v2xps_party ~ wth_prior + pi1*relevel(elections, ref = "Third-plus") + prior_pi + asp_dummy + + av_dm + (1 | ccode), data = wth8, REML=FALSE)

#get pi1 effect for each round of elections
tm1 <- tidy(m1, conf.int = TRUE) |> filter(term == "pi1") |> mutate(Sample = "None")

tm2 <- tidy(m2, conf.int = TRUE) |> filter(term == "pi1") |> mutate(Sample = "First")

tm3 <- tidy(m3, conf.int = TRUE) |> filter(term == "pi1") |> mutate(Sample = "Second")

tm4 <- tidy(m4, conf.int = TRUE) |> filter(term == "pi1") |> mutate(Sample = "Third-plus")

#now do it for prior_pi
m5 <- lmer(v2xps_party ~ wth_prior + prior_pi*elections + pi1 + asp_dummy + (1 | ccode), data = wth8, REML=FALSE)
m6 <- lmer(v2xps_party ~ wth_prior + prior_pi*relevel(elections, ref = "First") + pi1 + asp_dummy + (1 | ccode), data = wth8, REML=FALSE)
m7 <- lmer(v2xps_party ~ wth_prior + prior_pi*relevel(elections, ref = "Second") + pi1 + asp_dummy + (1 | ccode), data = wth8, REML=FALSE)
m8 <- lmer(v2xps_party ~ wth_prior + prior_pi*relevel(elections, ref = "Third-plus") + pi1 + asp_dummy + (1 | ccode), data = wth8, REML=FALSE)

#get pi1 effect for each round of elections
tm5 <- tidy(m5, conf.int = TRUE) |> filter(term == "prior_pi") |> mutate(Sample = "None")

tm6 <- tidy(m6, conf.int = TRUE) |> filter(term == "prior_pi") |> mutate(Sample = "First")

tm7 <- tidy(m7, conf.int = TRUE) |> filter(term == "prior_pi") |> mutate(Sample = "Second")

tm8 <- tidy(m8, conf.int = TRUE) |> filter(term == "prior_pi") |> mutate(Sample = "Third-plus")

#lastly we do it for asp
m9 <- lmer(v2xps_party ~ wth_prior + asp_dummy*elections + pi1 + prior_pi + (1 | ccode), data = wth8, REML=FALSE)
m10 <- lmer(v2xps_party ~ wth_prior + asp_dummy*relevel(elections, ref = "First") + pi1 + prior_pi + (1 | ccode), data = wth8, REML=FALSE)
m11 <- lmer(v2xps_party ~ wth_prior + asp_dummy*relevel(elections, ref = "Second") + pi1 + prior_pi + (1 | ccode), data = wth8, REML=FALSE)
m12 <- lmer(v2xps_party ~ wth_prior + asp_dummy*relevel(elections, ref = "Third-plus") + pi1 + prior_pi + (1 | ccode), data = wth8, REML=FALSE)

#get pi1 effect for each round of elections
tm9 <- tidy(m9, conf.int = TRUE) |> filter(term == "asp_dummy") |> mutate(Sample = "None")

tm10 <- tidy(m10, conf.int = TRUE) |> filter(term == "asp_dummy") |> mutate(Sample = "First")

tm11 <- tidy(m11, conf.int = TRUE) |> filter(term == "asp_dummy") |> mutate(Sample = "Second")

tm12 <- tidy(m12, conf.int = TRUE) |> filter(term == "asp_dummy") |> mutate(Sample = "Third-plus")

rbind(tm1, tm2, tm3, tm4, tm5, tm6, tm7, tm8, tm9, tm10, tm11, tm12) |> 
  mutate(term = case_when(term == "pi1" ~ "Incumbent PI", term == "prior_pi" ~ "Prior PI", term == "asp_dummy" ~ "ASP"),
         Sample = fct_relevel(Sample, c("None", "First", "Second", "Third-plus"))) |> 
  ggplot(aes(x = Sample, y = estimate, ymin = conf.low, ymax = conf.high, color = term, shape = term)) +
  geom_pointrange(position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(y = estimate - 100), linetype = "solid", color = NA) +  # Setting color to NA to remove legend
  geom_point(aes(y = estimate - 100), shape = NA, color = NA) +  # Setting color and shape to NA to remove legend
  theme_bw() +
  xlab("") +
  ylab("") +
  ylim(-0.55, 0.75) +
  scale_color_grey(name = "Variable", start = 0, end = 0.8) +
  scale_shape_manual("Variable", values = c(15, 16, 17)) + 
  theme(legend.position = "bottom")  +
  ggtitle("Random Effects Models") +
  theme(plot.title = element_text(hjust = 0.5))

#ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_fig1_2.jpg", dpi = 500, width = 5, height = 5)

#Table 1 of Appendix C----
wth2 <- wth |> select(v2xps_party, v2x_polyarchy, lgdp, wth_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()
wth3 <- wth |> select(countryname, year, v2xps_party, wth_prior, pi1, prior_pi, asp_dummy, ccode) |> na.omit()
wth4 <- wth |> select(v2xps_party, wth_prior, pi1, prior_pi, asp_dummy, wth_pduration, ccode) |> na.omit()
wth4 <- wth |> select(v2xps_party, wth_prior, pi1, prior_pi, asp_dummy, wth_pduration, regime_elecs, ccode) |> na.omit()
wth5 <- wth |> select(v2xps_party, wth_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, wth_number, fsu, ussr_sat, av_dm, ccode) |> na.omit()
wth6 <- wth |> select(v2xps_party, wth_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, wth_number, fsu, ussr_sat, av_dm, ccode) |> na.omit()
wth7 <- wth |> select(v2xps_party, wth_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, wth_pduration, wth_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()
wth8 <- wth |> select(ldv, wth_prior, pi1, prior_pi, asp_dummy, v2x_polyarchy, lgdp, wth_pduration, wth_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()
wth9 <- wth |> select(year, ccode, ldv, wth_prior, pi1, prior_pi, reactive, v2x_polyarchy, lgdp, wth_pduration, wth_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()

m1 <- lmer(v2xps_party ~ + (1 | ccode), data = wth, REML=FALSE) #

m2 <- lmer(v2xps_party ~ v2x_polyarchy + lgdp + wth_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = wth2, REML=FALSE) #

m3 <- lmer(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + (1 | ccode), data = wth3, REML=FALSE) # 

m4 <- lmer(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + wth_pduration + (1 | ccode), data = wth4, REML=FALSE) #

m5 <- lmer(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + wth_number + fsu + ussr_sat + av_dm + (1 | ccode), data = wth5, REML=FALSE)

m7 <- lmer(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + wth_pduration + v2x_polyarchy + lgdp + wth_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = wth7, REML=FALSE)

m8 <- lm(v2xps_party ~ wth_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + regime_elecs + av_dm, data = wth7) 
sd(residuals(m8))

m9 <- lm(ldv ~ wth_prior + pi1 + prior_pi + asp_dummy + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + regime_elecs + av_dm, data = wth8) 
sd(residuals(m9))

m10 <- lm(ldv ~ wth_prior + pi1 + prior_pi + reactive + v2x_polyarchy + lgdp + wth_pduration + wth_number + regime_elecs + av_dm, data = wth9) 
sd(residuals(m10))

modellist <- list("ML1" = m1, #just baseline
                  "ML2" = m2, #shows basline with time-variant indicators
                  "ML3" = m3, #just time invariant - prior duration 
                  "ML4" = m4, #add in prior duration, not much change 
                  "ML5" = m5, #full minus prior duration & # of elections, prior_pi drops 
                  "ML7" = m7,
                  "LM" = m8,
                  "LDV" = m9,
                  "LDV" = m10)

modelsummary(modellist, stars = TRUE, gof_omit = "R2|AIC|RM|F|Log.Lik", coef_map = c("pi1" = "Incumbent PI", "prior_pi" = "Prior PI", "asp_dummy" = "ASP",
                                                                                     "reactive" = "Reactive ASP",
                                                                                     "wth_priorPersonal" = "Personalist", "wth_priorMilitary" = "Military",
                                                                                     "wth_pduration" = "Prior Duration", "regime_elecs" = "Elections", 
                                                                                     "SD (Intercept ccode)" = "SD: Country", "SD (Observations)" = "SD: Residuals"), output = "latex")
#Figure 2.1 of Appendix C----
wth <- read_csv("data/wth_panel.csv") %>% 
  select(countryname, year, v2xps_party, v2x_polyarchy, wth_duration, wth_regime, wth_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, wth_pduration, wth_number, fsu, ussr_sat, regime_elecs, av_dm) %>% 
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown")) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(wth_prior %in% c("Personal", "Military", "One-party", "Monarchy") & is.nan(.) | wth_prior %in% c("Multi-party", "Military", "One-party", "Monarchy") & is.na(.), 0, .)) %>% 
  filter(wth_prior %in% c("One-party", "Military", "Multi-party") & wth_regime == "Democracy") %>% 
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) %>% 
  mutate(wth_prior = factor(wth_prior, levels = c("One-party", "Multi-party", "Military"))) %>% 
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) 

# Function to run the model with a specific baseline year
pi_models <- function(data, baseline_year) {
  # Relevel wth_duration to set baseline year
  data <- data %>% mutate(wth_duration = relevel(as.factor(wth_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm(v2xps_party ~ wth_prior + pi1 * wth_duration + prior_pi + asp_dummy, data = data)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> #filter(is.na(group)) |> 
    #mutate(ses = sqrt(diag(vcovCR(model, cluster = wth$ccode, type = "CR2"))), model = "None", conf.low = estimate - (qnorm(0.975)*ses), conf.high = estimate + (qnorm(0.975)*ses)) |> 
    filter(term == "pi1")
  
  return(coefficients)
}
prior_models <- function(data, baseline_year) {
  # Relevel wth_duration to set baseline year
  data <- data %>% mutate(wth_duration = relevel(as.factor(wth_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm(v2xps_party ~ wth_prior + prior_pi * wth_duration + pi1 + asp_dummy, data = data)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> 
    filter(term == "prior_pi")
  
  return(coefficients)
}
asp_models <- function(data, baseline_year) {
  # Relevel wth_duration to set baseline year
  data <- data %>% mutate(wth_duration = relevel(as.factor(wth_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm(v2xps_party ~ wth_prior + asp_dummy * wth_duration + pi1 + prior_pi, data = data)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> 
    filter(term == "asp_dummy")
  
  return(coefficients)
}

# Create an empty list to store the results
pi_results <- list()
prior_results <- list()
asp_results <- list()

# Loop through the first 10 years to use each as the baseline
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- pi_models(wth, baseline_year = as.character(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  pi_results[[as.character(year)]] <- coefficients
}
#prior models
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- prior_models(wth, baseline_year = as.character(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  prior_results[[as.character(year)]] <- coefficients
}
#asp models
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- asp_models(wth, baseline_year = as.character(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  asp_results[[as.character(year)]] <- coefficients
}

# Combine all results into a tidy dataframe
prior <- bind_rows(prior_results) 
pi <- bind_rows(pi_results) 
asp <- bind_rows(asp_results) 

bind_rows(prior, pi, asp) |> 
  mutate(term = case_when(term == "pi1" ~ "Incumbent PI", term == "prior_pi" ~ "Prior PI", term == "asp_dummy" ~ "ASP")) |> 
  ggplot(aes(x = as.factor(baseline_year), y = estimate, ymin = conf.low, ymax = conf.high, color = term, shape = term)) +
  geom_pointrange(show.legend = TRUE, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(y = estimate - 100), linetype = "solid", color = NA) +  # Setting color to NA to remove legend
  geom_point(aes(y = estimate - 100), shape = NA, color = NA) +  # Setting color and shape to NA to remove legend
  theme_bw() +
  xlab("") +
  ylab("") +
  ylim(-0.55, 0.75) +
  scale_color_grey(name = "Variable", start = 0, end = 0.8) +
  scale_shape_manual("Variable", values = c(15, 16, 17)) + 
  theme(legend.position = "bottom") +
  guides(color = guide_legend(order = 1), shape = guide_legend(order = 1))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_2_1.jpg", dpi = 500, width = 5, height = 5)

#Figure 2.2 of Appendix C-----
pi_models <- function(data, baseline_year) {
  # Relevel wth_duration to set baseline year
  data <- data %>% mutate(wth_duration = relevel(as.factor(wth_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm(v2xps_party ~ wth_prior + pi1 * wth_duration + prior_pi, data = data)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> #filter(is.na(group)) |> 
    #mutate(ses = sqrt(diag(vcovCR(model, cluster = wth$ccode, type = "CR2"))), model = "None", conf.low = estimate - (qnorm(0.975)*ses), conf.high = estimate + (qnorm(0.975)*ses)) |> 
    filter(term == "pi1")
  
  return(coefficients)
}
prior_models <- function(data, baseline_year) {
  # Relevel wth_duration to set baseline year
  data <- data %>% mutate(wth_duration = relevel(as.factor(wth_duration), ref = baseline_year))
  
  # Fit the model
  model <- lm(v2xps_party ~ wth_prior + prior_pi * wth_duration + pi1, data = data)
  
  # Extract the coefficient for pi1
  coefficients <- tidy(model, conf.int = TRUE) |> 
    filter(term == "prior_pi")
  
  return(coefficients)
}

# Create an empty list to store the results
pi_results <- list()
prior_results <- list()

# Loop through the first 10 years to use each as the baseline
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- pi_models(wth, baseline_year = as.character(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  pi_results[[as.character(year)]] <- coefficients
}
#prior models
for (year in 1:10) {
  # Run the model with the current year as the baseline
  coefficients <- prior_models(wth, baseline_year = as.character(year))
  
  # Add the baseline year to the result
  coefficients$baseline_year <- year
  
  # Store the result in the list
  prior_results[[as.character(year)]] <- coefficients
}

# Combine all results into a tidy dataframe
prior <- bind_rows(prior_results) 
pi <- bind_rows(pi_results) 

bind_rows(prior, pi) |> 
  mutate(term = case_when(term == "pi1" ~ "Incumbent PI", term == "prior_pi" ~ "Prior PI")) |> 
  ggplot(aes(x = as.factor(baseline_year), y = estimate, ymin = conf.low, ymax = conf.high, color = term, shape = term)) +
  geom_pointrange(show.legend = TRUE, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(y = estimate - 100), linetype = "solid", color = NA) +  # Setting color to NA to remove legend
  geom_point(aes(y = estimate - 100), shape = NA, color = NA) +  # Setting color and shape to NA to remove legend
  theme_bw() +
  xlab("") +
  ylab("") +
  ylim(-0.55, 0.75) +
  scale_color_grey(name = "Variable", start = 0.5, end = 0.8) +
  scale_shape_manual("Variable", values = c(16, 17)) + 
  theme(legend.position = "bottom") +
  guides(color = guide_legend(order = 1), shape = guide_legend(order = 1))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_2_2.jpg", dpi = 500, width = 5, height = 5)

#Figure 3 of Appendix C----
df <- read_csv("data/wth_panel.csv") |>  
  dplyr::select(countryname, year, v2xps_party, v2x_polyarchy, wth_duration, wth_regime, wth_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, wth_pduration, wth_number, fsu, ussr_sat, regime_elecs, av_dm) |>  
  mutate(ccode = countrycode(countryname, origin = "country.name", destination = "cown"), 
         reactive = case_when(ruling == 1 ~ 0, ruling == 0 ~ 1)) |> 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(wth_prior %in% c("Multi-party", "Military", "One-party", "Monarchy") & is.nan(.) | wth_prior %in% c("Multi-party", "Military", "One-party", "Monarchy") & is.na(.), 0, .)) |>  
  filter(wth_prior %in% c("One-party", "Military", "Multi-party") & wth_regime == "Democracy") |>  
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) |>  
  mutate(wth_prior = factor(wth_prior, levels = c("One-party", "Multi-party", "Military"))) |>  
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) 

#reactive*pi1##
library(interplot) #0.2.3
lm1 <- lm(v2xps_party ~ wth_prior + reactive*pi1 + prior_pi + v2x_polyarchy + lgdp + wth_pduration + wth_number + ussr_sat + regime_elecs + av_dm, data = df)

interplot(m = lm1, var1 = "reactive", var2 = "pi1", hist = TRUE, adjCI = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of a Reactive ASP on \n Party Institutionalization") +
  xlab("Incumbent PI") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold")) +  
  ylim(-0.9, 0.6)

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_3_1.jpg", dpi = 500, width = 5, height = 5)

library(marginaleffects) #0.25.1
lm2 <- lm(v2xps_party ~ wth_prior + reactive*pi1 + prior_pi + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + regime_elecs + av_dm, data = df)
me2 <- plot_slopes(lm2, variables = "pi1", condition = "reactive") 
dme2 <- me2$plot_env$dat

ggplot() +
  geom_pointrange(data = dme2, aes(x = reactive, y = estimate, ymin = conf.low, ymax = conf.high)) +
  ylim(-0.9, 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of Incumbent PI on \n Party Institutionalization") +
  xlab("Reactive ASP") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold"))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_3_2.jpg", dpi = 500, width = 5, height = 5)

#Figure 4 of Appendix C----
lm1 <- lm(v2xps_party ~ wth_prior + reactive*prior_pi + pi1 + v2x_polyarchy + lgdp + wth_pduration + wth_number + ussr_sat + regime_elecs + av_dm, data = df)

interplot(m = lm1, var1 = "reactive", var2 = "prior_pi", hist = TRUE, adjCI = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of a Reactive ASP on \n Party Institutionalization") +
  xlab("Prior PI") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold")) +  
  ylim(-1.2, 0.6)

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_4_1.jpg", dpi = 500, width = 5, height = 5)

lm2 <- lm(v2xps_party ~ wth_prior + reactive*prior_pi + pi1 + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + regime_elecs + av_dm, data = df)
me2 <- plot_slopes(lm2, variables = "prior_pi", condition = "reactive") 
dme2 <- me2$plot_env$dat

ggplot() +
  geom_pointrange(data = dme2, aes(x = reactive, y = estimate, ymin = conf.low, ymax = conf.high)) +
  ylim(-1.2, 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of Prior PI on \n Party Institutionalization") +
  xlab("Reactive ASP") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold"))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_4_2.jpg", dpi = 500, width = 5, height = 5)

#Figure 5 of Appendix C-----
lm1 <- lm(v2xps_party ~ wth_prior + reactive*pi4 + prior_pi + v2x_polyarchy + lgdp + wth_pduration + wth_number + ussr_sat + regime_elecs + av_dm, data = df)

interplot(m = lm1, var1 = "reactive", var2 = "pi4", hist = TRUE, adjCI = TRUE) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of a Reactive ASP on \n Party Institutionalization") +
  xlab("ASP PI") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold")) +  
  ylim(-0.3, 0.6)

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_5_1.jpg", dpi = 500, width = 5, height = 5)


lm2 <- lm(v2xps_party ~ wth_prior + reactive*pi4 + prior_pi + v2x_polyarchy + lgdp + wth_pduration + wth_number + fsu + ussr_sat + regime_elecs + av_dm, data = df)
me2 <- plot_slopes(lm2, variables = "pi4", condition = "reactive") 
dme2 <- me2$plot_env$dat

ggplot() +
  geom_pointrange(data = dme2, aes(x = reactive, y = estimate, ymin = conf.low, ymax = conf.high)) +
  ylim(-1.2, 0.8) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_bw() +
  ylab("Marginal Effect of ASP PI on \n Party Institutionalization") +
  xlab("Reactive ASP") +
  ggtitle("All Years Pooled") +
  theme(plot.title = element_text(hjust = 0.5), axis.title = element_text(size=10, face = "bold"))

ggsave(filename = "C:/Users/selfd/Documents/autopartyinstitutionalization/legacy/plots/appendix/appendix_c_5_2.jpg", dpi = 500, width = 5, height = 5)

#Table 2 of Appendix C------
wth <- read_csv("data/wth_panel.csv") %>% 
  mutate(COWcode = countrycode(countryname, origin = "country.name", destination = "cown")) |> 
  dplyr::select(countryname, year, COWcode, v2xps_party, v2x_polyarchy, wth_duration, wth_regime, wth_prior, prior_pi, pi1, pi2, pi3, pi4, asp_dummy, ruling, lgdp, wth_pduration, wth_number, fsu, ussr_sat, regime_elecs, av_dm) %>% 
  mutate_at(vars(pi1, pi2, pi3, pi4), ~ifelse(wth_prior %in% c("One-party", "Military", "Multi-party", "Monarchy") & is.nan(.) | wth_prior %in% c("One-party", "Military", "Multi-party", "Monarchy") & is.na(.), 0, .)) %>% 
  filter(wth_prior %in% c("Multi-party", "Military", "One-party") & wth_regime == "Democracy") %>% 
  mutate(asp_dummy = ifelse(is.na(asp_dummy), 0, asp_dummy)) %>% 
  mutate(wth_prior = factor(wth_prior, levels = c("Multi-party", "One-party", "Military"))) %>% 
  mutate(elections = ifelse(regime_elecs >= 3, "Third-plus", regime_elecs),  elections = recode(elections, `0` = "None", `1` = "First", `2` = "Second"), 
         elections = factor(elections, levels = c("None", "First", "Second", "Third-plus"))) |> 
  mutate(reactive = -1*ruling)

vdem <- vdem |> 
  dplyr::select(country_name, year, COWcode, v2xel_frefair)

df <- left_join(wth, vdem, by = c("COWcode", "year")) |> 
  rename(ccode = COWcode)

wth1 <- df |> dplyr::select(ccode, year, v2xel_frefair, v2xps_party, wth_prior, pi1, prior_pi, reactive, v2x_polyarchy, lgdp, wth_pduration, wth_number, fsu, ussr_sat, regime_elecs, av_dm, ccode) |> na.omit()

m1 <- lmer(v2xps_party ~ wth_prior + pi1 + prior_pi + reactive + wth_pduration + v2x_polyarchy + lgdp + wth_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = wth1, REML=FALSE)
m2 <- lm(v2xps_party ~ wth_prior + pi1 + prior_pi + reactive + wth_pduration + v2x_polyarchy + lgdp + wth_number + fsu + ussr_sat + regime_elecs + av_dm, data = wth1)
m3 <- lmer(v2xel_frefair ~ wth_prior + pi1 + prior_pi + reactive + wth_pduration + v2x_polyarchy + lgdp + wth_number + fsu + ussr_sat + regime_elecs + av_dm + (1 | ccode), data = wth1, REML=FALSE)
m4 <- lm(v2xel_frefair ~ wth_prior + pi1 + prior_pi + reactive + wth_pduration + v2x_polyarchy + lgdp + wth_number + fsu + ussr_sat + regime_elecs + av_dm, data = wth1)

sd(residuals(m2))
sd(residuals(m4))

modellist <- list("ML - PI" <- m1, #just baseline
                  "OLS - PI" <- m2,
                  "ML - Elections" <- m3, #just baseline
                  "OLS - Elections" <- m4) #shows basline with time-variant indicators)

modelsummary(modellist, stars = TRUE, gof_omit = "R2|AIC|RM|F|Log.Lik", coef_map = c("pi1" = "Incumbent PI", "prior_pi" = "Prior PI", "reactive" = "Reactive ASP",
                                                                                     "wth_priorOne-party" = "One-partyist", "wth_priorMilitary" = "Military",
                                                                                     "wth_pduration" = "Prior Duration", "regime_elecs" = "Elections", 
                                                                                     "SD (Intercept ccode)" = "SD: Country", "SD (Observations)" = "SD: Residuals"), output = "latex")
